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Abstract 

A novel approach to wall modeling for the incompressible Navier-Stokes equations including flows of moderate and 
large Reynolds numbers is presented. The basic idea is that a problem-tailored function space allows prediction 
of turbulent boundary layer gradients with very coarse meshes. The proposed function space consists of a standard 
polynomial function space plus an enrichment, which is constructed using Spalding’s law-of-the-wall. The enrichment 
function is not enforced but ’’allowed” in a consistent way and the overall methodology is much more general and also 
enables other enrichment functions. The proposed method is closely related to detached-eddy simulation as near-wall 
turbulence is modeled statistically and large eddies are resolved in the bulk flow. Interpreted in terms of a three- 
scale separation within the variational multiscale method, the standard scale resolves large eddies and the enrichment 
scale represents boundary layer turbulence in an averaged sense. The potential of the scheme is shown applying it to 
turbulent channel flow of friction Reynolds numbers from Re T = 590 and up to 5,000, flow over periodic constrictions 
at the Reynolds numbers Ren = 10,595 and 19,000 as well as backward-facing step flow at Re /, = 5,000, all with 
extremely coarse meshes. Excellent agreement with experimental and DNS data is observed with the first grid point 
located at up to y| = 500 and especially under adverse pressure gradients as well as in separated flows. 

Keywords: Wall modeling, LES, DES, turbulence, function enrichment, XFEM, variational multiscale method, 
law-of-the-wall, turbulent channel flow, periodic hill, backward-facing step 


1. Introduction 

Large-eddy simulation (LES) becomes prohibitively expensive for moderate and high Reynolds numbers if near¬ 
wall turbulence is resolved. Grid-resolution requirements enabling prediction of the necessary scales depend on the 
friction Reynolds number approximately as Re\ ID- Small computation cells in the boundary layer come along with 
severe constraints on the time step size to resolve the temporal scales of momentum-transfer mechanisms and to be 
compliant with the Courant-Friedrichs-Lewy (CFL) condition if explicit time integration schemes are utilized. 

The concept of wall modeling was therefore introduced in early works on LES of high-Reynolds-number flow by 
Deardorff El and Schumann 0 in an attempt to circumvent the resolution dependence on wall units. Wall modeling 
implies that near-wall turbulence and the accompanying momentum transfer are not resolved in detail but modeled in a 
statistical sense. With near-wall turbulence modeled, the size of dominating eddies in the bulk of the flow are governed 
by geometrical scales of boundary conditions with resolution requirements increasing approximately as Re 0A with the 
Reynolds number of the bulk flow Pfl . 

Common approaches in wall modeling focus on imposing synthetic boundary conditions prescribing tractions 
instead of no-slip velocity, see reviews e.g. in mm. This comes along with the major advantage that the velocity 
gradient does not have to be resolved explicitly by the scheme, enabling very coarse meshes. Yet accurate models are 
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required to predict the correct stresses. Simple equilibrium models allow for direct modeling of the tractions but are 
prone to inaccurate predictions in separated regions or flows with high pressure gradients ID. More accurate two-layer 
models have been developed, e.g. solving the simplified thin-boundary-layer equations (TBLE) on a separate domain 
between the wall and the first off-wall node to predict the momentum transfer inside the boundary layer 00121. For 
a comprehensive overview of wall-layer models it is referred to Piomelli and Balaras (4) and Piomelli |6). Alternative 
methods of imposing approximate boundary conditions have been proposed such as weak no-slip boundary conditions 
DMUl or the filtered-wall model m which are not frequently used but have lead to new insights in the field of wall 
modeling. 

Hybrid RANS/LES methods including detached-eddy simulation (DES) represent another paradigm for wall mod¬ 
eling mm. Instead of employing separate domains for near-wall and bulk turbulence such as in two-layer models, 
different subgrid models are applied in the respective regions of a single mesh. Reynolds-averaged Navier-Stokes 
equations are commonly employed in the wall region and LES subgrid closures in the bulk of the flow. The sharp ve¬ 
locity gradient present in high-Reynolds-number flow necessitates many computation points in wall-normal direction; 
under-resolution results in an imprecise prediction of the wall stresses leading to a log-layer mismatch 0. 

A major reason for DES requiring many grid points in wall-normal direction to be able to resolve the gradient 
accurately is the common application of low-order computational methods. In this work, we propose a problem- 
tailored high-order method in this region that is capable of resolving the velocity gradient with very coarse meshes. 
This is done by employing general knowledge about turbulent boundary layers without prescribing the velocity profile 
itself. Theoretical considerations on methods that allow for constructing customized numerical methods have first been 
introduced by Melenk and Babuska EO with their partition-of-unity method (PUM). Belytschko and Black |fl6l have 
subsequently suggested a formalism that allows for construction of a problem-tailored computational method in the 
application of crack propagation in solid mechanics. An enrichment function representing an approximate analytical 
solution is usually used to extend the solution space of the method, besides the standard polynomial function space. 
For a comprehensive overview of the method we refer to the review articles in El [El- 

Applications of this method in the field of fluid mechanics can be found in several academic examples such as 
enrichment with analytical high-gradient solutions of the convection-diffusion equation G2), simulating a sharp corner 
in Stokes flow via an asymptotic expansion as enrichment 1201 or resolving the bottom boundary layer of oceanic flow 
via a logarithmic enrichment function applied to a ID water column El- A recent publication suggests enrichment 
with modes obtained by a proper orthogonal decomposition to resolve the boundary layer of a stochastically forced 
Burger’s equation (22l . The general framework can also be used to resolve other features of the solution besides high 
gradients, such as jumps or kinks, and can even be used to cut elements (see, e.g., Ii23l ). In general, the enrichment 
function is not prescribed as solution but the method ’’chooses” the best solution among all functions available in its 
function space in a consistent manner. 

In this study, we suggest to extend the standard solution space with the law-of-the-wall due to Spalding as en¬ 
richment function such that the numerical method is able to represent the high-gradient velocity profile in a turbulent 
boundary layer with very coarse meshes. The idea follows the paradigm of DES as only the large eddies away from 
the wall are resolved explicitly and near-wall turbulent structures are represented in a statistical sense. The construc¬ 
tion of the method does however not require blending of turbulence models as the approach separates the statistical 
model a priori from the eddy-resolving space according to the variational multiscale method. Yet it is required that 
subgrid turbulence is modeled accurately such that the computational method is able to find an appropriate solution. 
Therefore a multifractal subgrid scale model embedded in a variational multiscale method is applied as subgrid scale 
model which gives excellent results for wall-resolved LES lf24l and has been successfully extended to passive-scalar 
mixing ||25| as well as low-Mach-number flow with variable density l26l . 

The present article is organized as follows. In Section [2] a weighted residual formulation of the incompressible 
Navier-Stokes equations is presented. Subsequently the enrichment of the function space with an appropriate gradient¬ 
capturing function is proposed in Section [3] and methods are presented for adaptivity in space and time. Turbulence 
modeling in the framework of the variational multiscale method and necessary adaptations to the new function space 
are revisited in Section [4] The method is validated with turbulent channel flow at moderate to moderately large 
Reynolds number, flow past periodic hills as well as flow over a backward-facing step in Section[5] Conclusions close 
the article in Section[6] 
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2. Incompressible Navier-Stokes equations 


The incompressible Navier-Stokes equations are considered in this work as outlined in the following. The 
weighted-residual formulation presented in the second subsection is the starting point for approximation of the solu¬ 
tion spaces in Section [3] 

2.1. Problem statement 

The incompressible Navier-Stokes equations in convective form are considered as 

— + it ■ Vn + Vp - 2vV ■ e(w) = / infix (0,T) (1) 

at 


V ■ u = 0 infix (0,7”) (2) 

with the fluid velocity u = (u\, U 2 , uf) T , the pressure p, the time t, the kinematic viscosity v and the symmetric rate- 
of-deformation-tensor e(u) = ^(V« + (Vu) T ). The body force vector is denoted /, the spatial domain Q and the 
simulation time T. At t = 0, a divergence-free initial velocity field is prescribed: 

u = uq in Pi x {0} (3) 

Dirichlet boundary conditions are defined as 

u-u D on T 0 x ((XT') (4) 

and traction boundary conditions are applied on the Neumann boundary 

crn-h onTjvXlO.T') (5) 

where the Cauchy-Stress tensor is cr = -pi + 2 ve(u). It is assumed that I f U Rv = 0 and \D Pi I ’v = T. 

2.2. Weighted residual formulation 

A weighted-residual formulation is obtained with a standard procedure by multiplying the momentum equation 
(|TJl with a weighting function v 6 *V y and the continuity equation (j2jl with q e 'V q . Appropriate spaces for u e S u . 
p e S p as well as *V v and 'V q are assumed. The choice of the discrete solution and weighting function spaces S„ 
and < V y is the main innovation presented in this article and is discussed in the subsequent Section [ 3 ] The equations 
are integrated over the domain Q, the pressure and viscous terms are integrated by parts and the Neumann boundary 
conditions are applied to the arising boundary integrals. The result reads 

S NS (v,q-,u,p) = l(v) (6) 

with the left hand side of the momentum equation and the contribution of the continuity equation 

du 

(V, q\ M, p) = (v, —) + (v, u ■ Vw) - (V • v, p) + (6(v), 2ve(ii)) + (q, V • u) (7) 

ot 

and the right hand side of the momentum equation 

l(v) = (v,f) + (v, h)r N - (8) 

The L2-inner product is defined as usual (•, •) = (•, -)n and (•, -)r v defines an integral over the Neumann boundary T lV . 
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3. Capturing the boundary layer via function enrichment 

The computational method applied to solve the incompressible Navier-Stokes equations has large influence on 
the quality of the solution and number of cells required. The method presented in the following utilizes a problem- 
tailored solution space distinguishing drastically from standard methods. The solution space of the method is capable 
of resolving high boundary layer gradients and adapts to local characteristics of the flow. This is done by extend¬ 
ing the solution space with the help of an approximate analytical representation of the mean velocity profile given 
as the law-of-the-wall. With a solution space capable of resolving the high gradient, the solution is not prescribed 
but the numerical method is able to find an appropriate solution in the offered function space, provided that turbu¬ 
lence is modeled accurately. The approach applied for turbulence modeling is discussed in Section [4] As numerical 
method, a variant of the PUM or the extended finite element method (XFEM) is chosen as it provides a framework for 
constructing such a customized space. 

3.1. Enriching the solution space 

The extended finite element method suggests a solution space u h (x, t) consisting of two contributions, a standard 
u h {x , t) and an enrichment part u h (x. t) as 


u h {x,t) = u h {x,f) + u\x,t) (9) 

dependent on the spatial location vector x — (xi , xi, X3) 7 and assuming direct sum decomposition of the underlying 
discrete solution spaces S„ = S„x S[\, which are identified with a characteristic element length h. The standard finite 
element expansion with shape functions and degrees of freedom u B is 

u h (x,t) = Y N u b (x)u b . (10) 

BeN u 

The enrichment including additional degrees of freedom m# is with the same partition of unity N 1 ^ defined as 

U h (x,t)= Y N b (x)(i//(x, t) - i//(x B , t))r'\x)u B (11) 

BeN“ tr 

where only a subset of nodes in the vicinity of the wall N" nr c N“ is enriched. The enrichment function t//(x, t) 
represents a problem-tailored profile, e.g. an analytical or approximate solution of the underlying problem and the 
enrichment function suggested in this article is presented in the following subsection. Subtracting the enrichment 
function by its nodal values i//(x B , t) yields zero on the nodes facilitating post processing and application of boundary 
conditions. 

Special treatment of the enrichment in the blending area is required. On the interface towards non-enriched 
elements, the enrichment does not vanish yielding Neumann boundary conditions with h = 0 on the enrichment 
nodes. As these nodes are subject both to in- and outflow of the domain and it is well-known that Neumann boundary 
conditions are ill-posed as inflow boundary for the incompressible Navier-Stokes equations, convergence problems 
are observed. This problem is circumvented by multiplying the enrichment with a ramp function r h (x) as detailed in 
Ii27l and depicted in figure [T] 

In a turbulent boundary layer, high-gradient solutions are only obtained for the velocity profile. The pressure space 
therefore remains unaltered and is given as a standard finite element expansion 

p h (x,t) = Y N p B (x)p B . (12) 

BeNp 

As partition of unity we choose shape functions of a standard eight-noded trilinearly interpolated hexahedral finite 
element throughout this article for N" and N p . Any other Lagrangian finite element could be employed as well 
including higher-order elements and unstructured grids. 
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Figure 1: Ramp function for blending. 



Figure 2: Decomposition of the mean velocity in a buffer layer into the linear and enrichment components. 


3.2. Law-of-the-wall as enrichment 

An appropriate choice of the enrichment function is the key feature of the overall methodology. This 

function provides the opportunity to include information a priori known about boundary layers in the function space 
without prescribing the solution itself. We propose to enrich the function space with an empirical single analytic func¬ 
tion for the law-of-the-wall including viscous sublayer and inner layer. This idea follows the paradigm of DES stating 
that not all turbulent scales need to be resolved at the wall but only their ensemble-averaged solution is computed. In 
contrast to standard methods for DES, the resolution in wall-normal direction may be very coarse if the function space 
is capable of resolving the mean gradient and the resolution requirements are essentially independent of wall units. 
The decomposition of the proposed function space into a linear standard and enrichment component is visualized in 
figure [2] 

Such mean velocity profiles have for example been suggested by Reichardt Ii28l or more widely known by Spalding 
f29l|. They satisfy the boundary conditions at the wall u (y = 0) = 0 and ^rr\ y =o = 1 (see e.g. Dean lf30l ) and may with 
the latter even predict the wall shear stress accurately. 

The enrichment function proposed in this work is a minor modification of the law-of-the-wall by Spalding 


y + (y, t w ) 


^ - 


<A 2 

1 — it/ — —— 
w 2 ! 


3! 4! 


(13) 


where the common formulation is recovered with u + = -. The constants k = 0.41 and B = 5.17 by Dean OH are 
applied. The only remaining parameters are the distance from the wall y and the wall shear stress t w included in the 
definition of the wall coordinate 


y + (y, Tw) 



(14) 


where the density is denoted p. Close agreement of Spalding’s law-of-the-wall with DNS data of turbulent channel 
flow at Re T = 2003 by Hoyas and Jimenez 021 is observed in figure[3]in the viscous sublayer and logarithmic region. 
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Figure 3: Comparison of the law-of-the-wall with DNS data of a turbulent channel flow at Re r = 2003. 


The modification of the law-of-the-wall by Dean lf30l to also describe the wake of the channel flow is not taken into 
consideration here. 

In the final method a discrete version y +,h (y h , rf y ) is used. The discrete distance from the wall y h is defined as the 
distance of each node ys, with B e N‘‘ nn to the closest node at the wall in /V" c N‘‘ nr 

y h = Yj n B y* (15) 

fle/V“ lr 

forming a robust procedure even for surfaces where the wall-normal vector is not unique. 

For application with the incompressible Navier-Stokes equations, expressions for the first and second derivatives 
of the enrichment with respect to cartesian coordinates are required. Their derivation is straight forward and included 
in|Appendix A| 


3.3. Adaptivity in space and time 

The wall shear stress being the only parameter of the function space represents an advantageous characteristic. It is 
known for many canonical flows in advance, such as for turbulent channel flow in a mean sense, and the pre-supposed 
value could be applied for computation of the wall model. However, in general flows, the tractions are not a priori 
known and it is desirable that the shape functions adapt to local fluctuations, regions of varying wall shear stresses 
and their temporal evolution. Therefore, it is suggested to explicitly compute the tractions and impose them on the 
wall model including spatial and temporal adaptation. 

Several methods are available for determining the instantaneous shear stresses, such as the common gradient-based 
approach via the derivative of the wall-parallel velocity component with respect to the wall distance, given in discrete 
form 


T 


h 

w 


£ »>|| 

BeN D 




(16) 


An alternative method common in the finite element method is calculating the wall shear stress via a nodal wall- 
parallel force vector r[j B on the Dirichlet boundary To divided by the nodally defined local area Ab and interpolated 
with the standard finite element expansion 


T 


h 

w 


R<=M U 


Ab 


(17) 


with the norm ||-[|2 of the three components corresponding to the space dimensions. The force vector equals the right- 
hand-side residual vector of the final matrix system and is discussed in Section 4.4 The nodal area is given as the 
integral of the standard partition of unity on the boundary 


Ab - 


f N“dA. 
J I'd 
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(18) 









Figure 4: Comparison of the wall shear stress (left) with the aggregated wall shear stress (right) of turbulent channel flow on a mesh with 32 3 
elements. Red indicates high and blue low values. 


Both of these methods represent an accurate definition of the momentary traction and are exactly equivalent for 
the continuous case but differences arise on discrete level. One of the differences is that the latter force-based method 
requires the residual to be converged to give an accurate prediction. Yet it is considered to be better consistent 
in the framework of FEM and chosen in this work as the standard method. The gradient-based method is applied for 
the first five time steps of the transient simulation as a converged residual is not available in the first time step and the 
gradient-based method is more robust if the initial field is not divergence-free. 

Another aspect that has to be considered in this context is that Spalding’s law is a relation for mean quantities, 
i.e. the mean velocity is related to the average wall shear stress. The difference between applying statistical and 
instantaneous values of the wall shear stress becomes apparent in the force-based method: The magnitude is computed 
for each node in equation ( fTT) resulting in a statistical over-prediction of the traction since neighboring force vectors 
usually are non-parallel. Therefore, it is suggested to calculate the stress via a locally averaged force field with a 
characteristic length scale ah instead of h resulting in a large-scale force. Such a local averaging operation allows for 
spatial variations of the traction and yet local fluctuations are smoothed. This averaging is realized via level-transfer 
operators from plain aggregation algebraic multigrid methods for separating scales, similar to the method frequently 
used to explicitly separate velocity scales in LES ||33| . A discrete wall shear stress r" h with a coarser characteristic 
element length ah as a multiple of the element length h is obtained. 

For this method a prolongation matrix P h ah is generated and the restriction matrix is defined as the transpose of the 
prolongation matrix resulting in R'f 1 = {P h ah ) T implying R" h h P l i l rh — / with the identity matrix /. A scale-separation or 
aggregation operator is defined as 

S ? = ph ah R T (19) 

yielding a coarse-scale force field via a vector-matrix multiplication 


>,ah 


_ r <ah v 

~ r 


( 20 ) 

This result is applied to calculate the shear stress r" v h in equation (17 1 with a value of a = 3. Figure [ 4 ] compares r* 
and showing that the field variable is averaged locally for the latter one but still may take larger variations into 
account. Thus, is an appropriate representation of the wall-shear stress for spatial adaptation of Spalding’s law. 

The traction computed at the no-slip nodes is communicated to the respective off-wall nodes in N l ‘ m . where the 
node pairs are again determined via the shortest distance. 

The temporal evolution of the turbulent flow results in continuous adaptation of the function space. For simplicity, 
the space is kept constant during the non-linear Newton iterations implying a quasi-static treatment. Updating the wall 
shear stress yielding new shape functions requires a subsequent L2-projection of the solution of the previous time step 
n onto the new space of the current time step n + 1 given in weak form as 


(-A.n+1,-/!,»+!) = ^h,n+\ a h,ny 


( 21 ) 


Complementary vectors required by the discrete time integration procedure such as a potential acceleration vector are 
projected with the same matrix system. 
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4. Subgrid-scale modeling 

The major goal of all wall models is to utilize very coarse meshes in the near-wall region implying that a large part 
of the physics is not resolved but modeled. This defines a need for accurate turbulence models ensuring high-quality 
results. In contrast to the common approach in LES of using a filter in order to derive the unresolved part ||34& , a scale 
separation as suggested originally for large-eddy simulation by Hughes et al. |[35l is used here. The scale separation 
gives rise to unresolved scales that we model by a structural reconstruction via a multifractal subgrid scales turbulence 
model embedded in a residual-based variational multiscale method as outlined in the following. 

4.1. Scale separation for large-eddy simulation 

In addition to the decomposition of the solution space into standard and enriched components, the velocity space 
is further separated into resolved u h and unresolved scales ii reading 

u = u h +u- u h +u h + u. (22) 

The resolved scales are again identified by a characteristic element length h and the unresolved scales with (•). The 
equivalent separation of scales of the pressure into a resolved p h as well as unresolved component p is also performed 
resulting in 

P = P h + P- (23) 

Direct sum decomposition of the underlying resolved and unresolved spaces is assumed as S u = S„ xS u = xS„ xS„ 
and S p - Sp x S p . Inserting these definitions into the weighted residual formulation (jtijl gives rise to the following 
terms: 

& NS (v,q;u h ,p h )+B l f' s (v,q;u,p) + C(v;u l \u) + <R(v;u) = Z(v) (24) 

The term &Ns(v,q\ u h ,p h ) constitutes the part of the relation that is represented by the resolved solution space. The 
contribution B l ™ s (v, q\u,p) summarizes the linear terms dependent on the subgrid scales u and p: 

S'f' s (v,q-u,p) = (v, ^) - (V ■ v,p) + (e(v),2 ve(u)) + (q,V ■ u) (25) 

The non-linear convective term gives rise to the cross- and Reynolds stresses as 

C(v; u'\ u) = (v, u h • Vm + u ■ Vm /! ) (26) 

and 

7?(v; u) - (v.u ■ Vm). (27) 

A basic characteristic of the variational multiscale method is that the solution and weighting function spaces have 
the same structure, i.e. the spaces for the weighting function may likewise be decomposed into the corresponding 
resolved and unresolved contributions: 

V = V h + v = v h + v h + V (28) 




q h +q 


(29) 


Since equation (24 1 is linear with respect to the weighting functions, it may be separated and as usual only the resolved 
scale component is taken into further consideration: 


S NS (v l ',q h -u l ’,p h ) + £% (v\q h ;u,p) + C(v h ;u h ,u) + <R(v h ;u) = l(v h ) 


(30) 


This result still contains the unresolved scale quantities u and p which are unknown and have to be modeled. In 
the following sections, this is done via multifractal scale similarity as well as residual-based modeling. 

Remark: The enrichment approach presented in this paper may also be interpreted as a separation of the solution 
vector in three scale groups as for example described by Gravemeier et al. 
three scales are represented by the standard resolved scale u h 


and indicated in relation (22 1 . The 
the enrichment scale u h as well as the unresolved scale 


h. With regard to LES, the standard scale resolves large eddies that are at least of the size of the characteristic element 
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length. The enrichment scale, however, represents flow features in a statistical sense and without resolving large eddies 
in the near-wall region explicitly. The physical interpretation of the unresolved scales are fluctuations on subgrid level. 
This interpretation and the framework of the variational multiscale method would allow different turbulence models 
tailored for the respective requirements of the scale. The turbulence models presented in the following show to be 
sufficiently general and adapt to the specific requirements in different regions of the domain such that a three-scale 
turbulence model is not necessary. 


4.2. Subgrid modeling with multifractal subgrid scales 

The cross and Reynolds stress terms ( [26] ) and ( [27] ) are modeled explicitly by reconstruction of the unresolved scale 
it via a multifractal subgrid-scale model as proposed by Rasthofer and Gravemeier 11241 . The multifractal subgrid scale 
model follows the idea that turbulence originates from repeated stretching and folding of vortical structures and that 
this process is scale-invariant. The model attempts to reconstruct the subgrid-scale vorticity and computes the subgrid 
velocity through the law of Biot-Savart, indicating that the large eddies of the flow have to be resolved explicitly. For 
a detailed derivation of the governing relations it is referred to Burton and Dahm (37). 

The subgrid velocity scales with the small-resolved velocity field 6u h and a proportionality factor B as 

it ~ B6u h . (31) 


The small-scale velocity is determined by an explicit filtering of the standard FE part of the resolved velocity yielding 


u = u ah 


■ 6u h + u h + u 


(32) 

for the overall composition of the velocity space. The large-scale velocity field u ah is identified by a length scale of 
ah as a multiple of the element length. This decomposition is chosen due to the physical interpretation of the standard 
finite element space as eddies while the enrichment space represents a statistical velocity profile that does not resolve 
eddies by nature. 

Explicit scale separation of u ah and 5u h is performed via a plain aggregation algebraic multigrid method as pro¬ 
posed in j33l and applies similar relations as used for smoothing of the wall shear stress in Section 3.3 The standard 
parameter a = 3 is applied. 


The proportionality factor in (31 1 is given as 


4 1 2N 4N 

B = C sgJ (l - a 5) 22 3 (2 3 


1)5 


(33) 


In the current application of convection-dominated high-Reynolds-number flow, the constant C sgs = 0.15 is chosen. 
This value is significantly lower than the one suggested in l24l . which is since the near-wall limit as suggested in the 
original publication is not considered here. B is evaluated at the quadrature points during evaluation of the discrete 
formulation ( |30] >. The number of cascade steps N from the smallest resolved scales of size h to the viscous scale A,, is 
approximated via the local element Reynolds number Ret, = ^ lh and a proportionality constant c v resulting in 


N = log 2 (—) = log 2 (,c v Re i h ). 
A v 


(34) 


A value for the proportionality constant c v = 0.1 is used, which is close to the value of determined experimentally 
by Mullin and Dahm l38l and h is approximated by the cube root of the local element volume. 


The final result for the modeled cross and Reynolds stresses (26 1 and (27 1 with the presented relation for the 
subgrid-scale velocity © is 

C(v h - u\ u) ~ (v h , B(u h ■ V6u h + 6u h ■ V« A )) (35) 


and 


7?(v"; u) 


(v h , B\6u h ■ V6u h )). 


(36) 
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4.3. Residual-based modeling 


The multifractal subgrid scale model as presented in the previous section enables reconstruction of the subgrid 
velocity field and by that models the cross and Reynolds stress terms. As reported in (39), the model allows both for 
dissipation and backscatter of energy resulting in potentially de-stabilizing effects. Therefore, as suggested in i24l.it 
is embedded in the residual-based variational multiscale method providing a stable numerical method. 

The remaining linear terms of the scale separation (25 1 are approximated with 


&™(y h ,q h -,u,p) * (u h V • v* T M R h M ) + (V • v h ,T C R h c ) + (V, T M R h M ) (37) 

SUPG grad-div PSPG 


The included terms consist in a Streamline/Upwind Petrov-Galerkin (SUPG), a grad-div (grad-div) and a Pressure 
Stabilizing Petrov-Galerkin (PSPG) term. The SUPG term stabilizes the method with respect to convection by in¬ 
troducing a certain amount of artificial dissipation BUI . Better fulfillment of the divergence-free constraint (|2]i and 
improved convergence of the iterative solver is obtained via the grad-div term ED which also introduces a certain 
amount of dissipation in the system. The PSPG contribution enables circumventing the inf-sup condition (see e.g. 
J42|) and allows equal-order interpolation J43j. 

The momentum residual R f ‘ is defined as 


du h 

dt 


R% = — + u" ■ Vh" + Vp" - 2vV • e(u ') - f + B(u h ■ Vdu” + 6u" ■ Vh") + B~(6u • V<5h"). (38) 


In contrast to [24], it is suggested to include the modeled cross- and Reynolds stress terms (261 and (27 1 in the residual 
for better consistency. The discrete continuity residual R 1 ' is 


= V • u h . 


(39) 


The stabilization parameters tm and t c are designed to take into account the non-polynomial character of the 
element space. A definition inspired, among others, by Codina 04) and Gravemeier et al. j45l is chosen including, as 
usual, a transient, convective and viscous contribution for tm as 


1 

t m - - i= - 

l t+ 2 y ]f\\u%+4A l, y 


(40) 


with the time step At and a reciprocal scaling of t m and r c yielding 


1 


TC 


4 A h Th 


(41) 


It is noted that ^2 yj -< 4 A h which has been reported to be a requirement for example in |[44 ]. 

The parameter A h generally incorporates the characteristics of the element, for example the polynomial order of 
the underlying function space. For standard Lagrangian elements with well-defined polynomial order, such as the 
non-enriched elements, values are for instance for the polynomial orders p = {1,2,3} given as A h = {^, } with 

the characteristic element length h 14611471 . In the current application, element spaces of the enriched elements are 
non-polynomial making it impossible to specify appropriate values a priori. Therefore, an element-specific value 
for A h is determined consistently via inverse estimate as suggested by Harari and Hughes |[46l ensuring stability for 
convection-dominated flows by solving the local generalized eigenvalue problem for its maximum eigenvalue A h and 
v h given as 


(AwAAvV - AVwAVvV = 0. 


(42) 


Q'' represents the element domain and the solution v h and weighting function space \v h are defined similarly as the 
enriched velocity space (|9|, e.g. for v h : 

v h (xj) = v h (x,t) + v h (x,t) (43) 
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The standard and enrichment parts are given with a single degree of freedom per node as 


A x,t) = ^ Ng(x)v B 


(44) 


and 


v h (x, t) = ^ N“(x)(tJ/(x, t) - \p(x B , t))r h (x)v B . 


(45) 


BeN“ 


Evaluating the inverse estimate with only one degree of freedom per node results in matrix dimensions of only 
16 x 16 for the elements presently considered, such that the eigenvalue-related computation time is negligible. A great 
characteristic of the presented stabilization parameter is highlighted: tm and t c are completely free of the element 
length if A h is determined via (42 1 . Especially for anisotropic elements, the definition of h is not obvious and many 


definitions have been proposed. The advantages of such a definition are also discussed for example by Franca and 
Madureira j48l . 

Solely for linear elements, the standard value of A h = ^ is applied with the volume-equivalent diameter h = 
(—)s / s/PT^j with V the element volume and the number of space dimensions n v ,/ for simplicity (491. Due to large 
variations of the stabilization parameters tm and Tc within one element, especially in the first element at the no-slip 
boundary condition, the parameters are evaluated at the quadrature points. 

The present approach for turbulence modeling has been presented as a subgrid-scale model for LES assuming 
that the largest eddies present are resolved by the scheme. However, inside the first element at the wall, with the first 
off-wall node placed at y + > 100, this is certainly not fulfilled in the wall-region. The presented turbulence model has 
yet proven to be able to model the necessary subgrid scales even if the largest eddies are not resolved everywhere. 


4.4. Final discrete problem 

The final semi-discretized problem becomes 

bu h 

(A —) + (v\ u h • Vm ,! ) - (V • v\ p h ) + (e(v' ! ), 2ve(u h )) 
ot 

+ ( v h , B(u h ■ V 6u h + 8u h ■ Vu h )) + (A B 2 (6u h ■ V<5m /! )) 

C K 

+ (u h V-v h ,T M R h M ) + (V • v\ T (: R b c ) 

SUPG grad-div 

+ (Av-iA + (vAtm<) 

PSPG 

= (A/*) + (A h h \ N (46) 

where the contributions of multifractal and residual-based subgrid-scales are labeled. The residuals R h M and R h ( . are 
defined in ( [38] ) and ( |39| ), the stabilization parameters t m and t c are given in ( |40] ) as well as ( |4~T[ i and B in ( |33[ >. The 
terms are integrated in space applying direction-dependent GauB-quadrature rules of appropriate order that enable 
accurate integration despite the non-polynomial function space. Equation ( |46[ is integrated in time utilizing a second- 
order accurate generalized-rr time integration scheme including= 0.5 1150115II . Adaptive time stepping is employed 
such that the maximum CFL condition is kept constant at CFL=0.5 for all simulations presented. 

The final matrix system is linearized and iteratively solved via a Newton-Raphson scheme yielding 


OT+li J+l _ H +1 
ZAZ i+1 - r i 


(47) 


for the current time step n + 1 and non-linear iteration i + 1, omitting the superscript h for simplicity. The increment 
includes both velocity and pressure increments from the current non-linear iteration such that 


U i+ 1 - Vi 

Pi + 1 - Pi 


Az/+i = 


A£/,+i 

AP i+1 
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Case 

N x i x N x2 x N x3 

iaoie i: L.nannei now cases ana resolutions. 

Re T 

y* 

N wm 

Ch8wm2 

8 x8x8 

590; 950; 2,000 

147.5; 237.5; 500 

2 

Ch\2wm2 

12 x 12 x 12 

590; 950; 2,000 

98.3; 158.3; 333.3 

2 

Chl6wtn3 

16 x 16 x 16 

590; 950 

73.8; 118.8 

3 

Ch\6wm2 

16 x 16 x 16 

2,000 

250 

2 

Ch2Awnii 

24 X 24 X 24 

950 

79.2 

3 

Ch2Awm2 

24 x 24 x 24 

5,000 

416.7 

2 

Ch32wm2 

32 x 32 x 32 

5,000 

312.5 

2 

Ch96 

96 x 96 x 96 

950 

1.4 

- 



Figure 5: Decomposition of the mean velocity u + = u\/u T with u T = of case Ch\1wm2 at Re T = 2,000 into linear and enrichment part and 

comparison to DNS data. Symbols indicate nodes and the first off-wall node is located at y| = 333.3. 


The matrix K contains the linearization of all contributions of ( |46| ) except C, R and the respective stabilization terms, 
which are treated in a fixed-point-like procedure 1241 . The residual r summarizes all terms of ( [46] ) at the previous 
non-linear iteration i. In ( |49] > K is split into four parts including K' u , K vp , K qu and K qp and r is split into two vectors 
r v and r q : 

\K VU K vp 

^ ~ k>i" [(‘ip 



K vu contains the transient, convective and viscous term as well as terms of SUPG and grad-div. K vp comprises the 
pressure term and the respective part of SUPG. K qu and K qp summarize the continuity contribution and the PSPG 
terms. The nodal values of the momentum-residual vector r v on the Dirichlet boundary are equivalent to the nodal 
forces and are used to calculate the wall-shear stress r" h in equation (17 1 . 


5. Numerical examples 

In this section, the performance of the presented approach is investigated for turbulent channel flow at various 
Reynolds numbers, flow over periodic hills and backward-facing-step flow. The latter two examples discuss the 
performance under separated boundary layer conditions and adverse pressure gradients. 

5.1. Turbulent channel flow at moderate and moderately large Reynolds numbers 

A channel of the dimensions 2nd x 2d x nd in streamwise, wall-normal and spanwise direction, respectively, with 
periodic boundary conditions and channel-half height 6 is considered. We discuss flows of friction Reynolds numbers 
Re T = 590, 950, 2,000 and 5,000 on very coarse uniform meshes with 8 x 8 x 8 up to 32 x 32 x 32 elements, see 
table[T]for an overview. To the best of the authors’ knowledge, turbulent channel flow of friction Reynolds numbers 
higher than Re T — 950 has so far never been published employing a residual-based turbulence modeling approach. 
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Figure 6: Normalized mean velocity for Re T = 590, 950, 2,000 and 5,000, each shifted upward by 6 units for clarity. Symbols indicate nodes. 


The column N wm indicates the number of layers of enriched elements employed at the solid boundaries. Table[l]also 
compares the location of the first off-wall nodes in wall units y|, which are located between yf = 73.8 and yf = 500 
for the different flows and given discretizations. For comparison, we also include a simulation with resolved near-wall 
region without wall model at Re T = 950 on a discretization with 96 x 96 x 96 elements. The results presented in the 
following are labeled according to table [T] They are compared to direct numerical simulation (DNS) data of Moser, 
Kim and Mansour l52ll for Re T = 590, Del Alamo and Jimenez ll53l for Re T = 950 and Hoyas and Jimenez lf32l for 
Re r = 2,000. The results for Re T = 5,000 are compared to u + = -ln(y + ) + B with k and B defined as in Section ’ 


3.2 


We commence the discussion of the results with figure [5] showing the decomposition of the mean velocity of case 
ChYlwml at Re T = 2,000 similar to figure j^J The normalized mean velocity profile u' = ^ with u T = follows 
DNS data closely and provides an excellent match despite the extremely coarse resolution. With the first off-wall node 
located at y| = 333.3, a large part of u + is in the first element represented by the enrichment part of the flow it' which 
also constitutes the largest part of the gradient at the boundary. Away from the wall in the second element layer, the 
contribution of the standard space u + constitutes almost the whole solution. 

The non-dimensional mean velocity for all 13 simulations included in table [T] is shown in figure [6] A striking 
independence of the mesh applied is observed for all Reynolds numbers. Even for discretizations consisting of only 
8x8x8 elements, where hardly the largest eddies are resolved, the results are in acceptable to very good agreement 
with DNS data. Further, it is noticed that the normalized mean velocity is slightly over estimated for friction Reynolds 
numbers Re T = 590 and 950. The simulation results at Re T = 2,000 and 5,000 exhibit an excellent match with 
reference data, however. Comparing the wall-resolved LES of case Ch96 at Re r = 950 a mean velocity of slightly 
lower quality than the wall-modeled LES is obtained. Very good agreement with LES data is first obtained with a 
mesh of 128 x 128 x 128 elements as shown by Rasthofer and Gravemeier 11241 . 

The performance of the present wall model is further assed via root-mean-square (RMS) data of the fluctuations 
m' + , v’ + in wall-normal as well as w' + in spanwise direction and Reynolds shear stresses (t<V) + at Re T = 950 displayed 
in figure [7] Considering u ,+ , a distinct tendency to convergence for an increasing number of elements is observed. 
For w' + , the predictions show a similar behavior as observed for u’ + while v' + is generally predicted too small. The 
Reynolds shear stresses are predicted quite accurately for all discretizations. That near-wall fluctuations are not of 
the same quality as mean velocities is presumably inherent in the wall-modeling approach as the enrichment function 
constitutes a mean-velocity profile and the major part of the fluctuations is not resolved due to the coarse meshes 
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Figure 7: RMS of velocity fluctuations u' + = rms(u\ )/u T , v' + = rms{u' 2 )/u T and w' + = rms(u\)j u T as well as Reynolds shear stresses (u'v') + = 
(«l «/)/«/ for Re T = 950. 


applied. As expected, RMS data of the case C7t96 is in favorable agreement with DNS data as a significant amount of 
near-wall turbulent structures is resolved. 

From the results presented in this section it is concluded that extremely coarse resolutions can be used for simula¬ 
tion of turbulent channel flow. The first off-wall node can be located at up to yf = 500 and striking results are obtained 
even for meshes consisting of only 8x8x8 where hardly the largest eddies are resolved. 

5.2. High-Reynolds-number flow over periodic constrictions 

We consider flow over a smoothly curved 2D-periodic hill as described and analyzed numerically e.g. by Frohlich 
et al. f54l and Breuer et al. l55l and experimentally by Rapp Il56l with a Reynolds number based on the hill height H of 
Ren = 10,595 as well as 19,000 to validate our wall modeling approach. This flow configuration includes separation 
at the crest of the hill, a recirculation bubble formation, reattachment as well as recovery. Wall models for LES are 
challenged by this flow with a strong adverse pressure gradient that causes many models to produce deficient results. 
For example, Chen et al. 0 have found that their wall model based on the simplified TBLE under estimates the skin 
friction in the recirculation region as the convective term is neglected in that model. In the current wall modeling 
approach, all terms of the Navier-Stokes equations are retained such that a better performance with respect to adverse 
pressure gradients may be expected. Temmerman et al. ED investigated several wall functions and subgrid closures 
for LES and found that the location of the separation point has major impact on the reattachment location. Also, 
accurate prediction of the separation point of this flow are challenging employing steady RANS simulations ||58l . 
Hybrid RANS/LES techniques have been analyzed by Breuer et al. Il59l and Saric et al. l60l who have shown that 
the RANS/LES interface should be located inside the boundary layer on the crest of the hill. Due to the construction 
of the present DES technique, there is no explicit interface between the statistical and the LES region such that these 
problems are not expected to occur. 

A domain of the dimensions 9 H x 3.036 H x 4.5 H with periodic boundary conditions in the streamwise and 
spanwise direction and no-slip boundary conditions at top and bottom is considered. A very coarse mesh comprising 
64 x 32 x 32 as well as a refined, yet very coarse, mesh with 96 x 48 x 48 cells with uniform grid spacings in all 
directions and vertical grid lines as depicted in figure[8]are utilized. An overview of the simulations presented is given 
in table [2] where the coarser grid is labeled PhC and the finer grid PliF for Re H = 10,595. A simulation without 
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Table 2: Simulation cases and resolutions of the periodic hill. Ren = 10,595: PhC coarse mesh with wall modeling, PhF refined mesh with wall 
modeling, PhFNWM refined mesh without wall modeling, Froehlich.et.al. highly resolved LES, RappJEXP experiments, Chen_et_al._C and 
Chen_et.al.-F wall modeling based on simplified TBLE and immersed interface method. Ren = 19,000: PhC 19 coarse mesh with wall modeling, 
PhF\9 refined mesh with wall modeling, PhFNWM\9 refined mesh without wall modeling, RappJEXP 19 experiments. 


Case 

N x \ X N x 2 x N x3 

Re H 

X\,sep! H 

X\ ,reattl H 

Nwm 

PhFNWM 

96 x 48 x 48 

10,595 

0.2 

3.68 

- 

PhC 

64 x 32 x 32 

10,595 

0.25 

3.77 

4 

PhF 

96 x 48 x 48 

10,595 

0.25 

4.91 

4 

Froehlich_etMl. 

192 x 128 x 186 

10,595 

0.2 

4.6 - 4.7 

- 

Rapp-EXP 

- 

10,600 


4.21 

- 

Chen_et_al.JC 

96 x 64 x 32 

10,595 

0.65 

4.0 

- 

Chen etjal.-F 

192 x 72 x 48 

10,595 

0.5 

4.42 

- 

PhF NWM 19 

96 x 48 x 48 

19,000 

0.2 

3.4 

- 

PhC\9 

64 x 32 x 32 

19,000 

0.24 

2.58 

4 

PhF 19 

96 x 48 x 48 

19,000 

0.26 

3.94 

4 

Rapp EXP\9 

- 

19,000 


3.94 

- 




Figure 8: Grid of case PhC. Enriched elements are colored red and 
standard linear elements blue. 


Figure 9: Location of the first off-wall node in wall units over the 
streamwise coordinate of the periodic hill case. 


wall model is also investigated for comparison, which is labeled PhFNWM and employs the finer mesh. Considering 
Ren = 19,000, the same meshes are applied and labeled PhC 19, PhF 19 and PhFNWMl9, respectively. The number 
of enriched element layers is N wm = 4 for both meshes on both top and bottom wall. Figure [9] shows the location of 

the first off-wall grid point estimated as y|‘ / ' = y over the X| -coordinate. Here, y / ' is not the actual wall distance 
but the distance to the closest node at the wall. The first off-wall node is located at varying distance depending on 
resolution and Reynolds number up to approximately y^’ h =216 with minima near the zero-crossings of the wall shear 
stress. The mass flow is kept approximately constant over the simulation time and statistics are sampled over 10,000 
time steps. 

An overview with respect to reference data considered is also given in table [2] The results for Re H - 10,595 are 
compared to data of highly resolved LES by Frohlich et al. ED, labeled as FroehlichMjal. as well as the coarse mesh 
discussed by Chen et al. 0 (>Chen.etMl.JC ) with 64 cells in vertical direction. The separation and reattachment points 
are further compared to experiments by Rapp |56l ( Rapp_EXP ) and the fine mesh by Chen et al. 0 ( Chen_et_al.-F ). 
The results for Ren = 19,000 are compared to experiments by Rapp ll56ll (Rapp_EXP 19). Data labeled as Rapp_EXP 
and Rapp_EXP\9 has been obtained from the ERCOFTAC QNET-CFD Wiki contributed by Rapp et al. |61|. 

We start with a discussion of the results for the flow at Ren = 10,595. The skin-friction coefficient c/ at the lower 
wall and pressure coefficients c p at the upper and lower wall are compared to resolved LES data in figure [TO] They 
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Figure 10: Skin-friction (left) and pressure (right) coefficients for the flow over periodic constrictions at Ren = 10,595. The shallower pressure- 
coefficient curves correspond to the top wall. 


are defined as 


and 


II 

0 

(50) 

iP lt b 


P - Pref 

C P ~ 1 ,..2 

(51) 


iP u b 


with the bulk velocity if, and r w via the right-hand side residual (17 1 . As reference pressure p re f , the pressure at the 
upper wall at x\ =0 is chosen. The skin friction computed via wall-modeled LES is in close agreement with reference 
data over large parts of the domain. Solely at the crest of the hill, the peaks between x\ = 8 and 9 as well as x\ =0 and 
1 are significantly over-predicted, but improve for the case PliF with higher resolution. This over-prediction might 
be related to the local averaging operation of the wall shear stress applied during construction of the shape functions. 
The minor recirculation at the top of the hill observed in highly resolved LES data is not visible in the results of PhC 
and PhF due to the coarseness of the mesh. Separation and re-attachment points are predicted accurately via the zero¬ 
crossing of cj as well and are summarized in table[2] For the case PhFNWM without enrichment, large discrepancies 
including high peaks are visible on top of the hill. In the recirculation region, the skin-friction is over estimated and 
the reattachment length is predicted shorter than for the cases with wall modeling. The skin-friction coefficient is also 
compared to the results of Chen et al. a who under estimate Cf significantly due to the neglected convective term as 
aforementioned. 

The pressure curves of the present wall model are also in very good agreement to reference data and improve 
with resolution. The case PhC shows minor discrepancies both at the lower and upper wall, which are due to the 
coarseness of the resolution. The case PhFNWM over-predicts the pressure in the recirculation bubble and exhibits 
negative peaks on the hill crest. In the recovery region, the estimation is comparable to the coarse mesh with wall 
modeling PhC. At the upper wall, the prediction with wall model is superior compared to the one without. 

Mean velocities u\ in streamwise and 112 in vertical direction and Reynolds-shear stresses u\ u' 2 and the turbulent 
kinetic energy (TKE) k = i(n' 2 + u 2 + it'?) of the case Re H = 10,595 at ten stations are compared with the LES data 
of Frohlich et al. Il54l in figure [IT] The mean velocity u\ exhibits discrepancies with reference data for case PhC in 
the reattachment and recovery region while the finer mesh PhF results in a perfect match with reference data. Without 
wall model, PhFNWM predicts u\ with similar quality as the coarse mesh with wall modeling, PhC. Also for the 
mean velocity 112 in vertical direction, excellent results are obtained for the fine mesh including wall modeling PhF. 
For the other cases, M 2 is under estimated above the recirculation bubble due to the shorter re-attachment length. The 
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PhFNWM - PhC 


PhF 



ui/ut + xi/H 





Figure 11: Mean velocity u\ in xi and 1/2 in at- direction, turbulent kinetic energy k as well as Reynolds shear stresses u' u' 2 for the periodic hill at 
Refj = 10 , 595 . 
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Figure 12: Mean velocity »i in xi and u 2 in Ay-direction as well as Reynolds shear stresses «' u' 2 for the periodic hill at Ren = 19,000. 
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Table 3: Simulation cases and resolutions of backward-facing step flow. BFS-NWM without wall modeling, BFS-WM 3 with wall modeling, 
BFS J&.D.EXP experiments, BFS LMK DNS DNS. 


Case 

Re h 

X\,reatt! h 

N wm 

BFS _NWM 

5,000 

13.49 

- 

BFS JVM3 

5,000 

6.78 

3 

BFS _J&D_EXP 

5,000 

6.0 + 0.15 

- 

BFS LMKJDNS 

5,100 

6.28 

- 




Figure 13: (Top) Instantaneous velocity magnitude over the backward-facing step: red indicates high and blue low values. (Bottom) Mesh in the 
vicinity of the step: enriched elements are colored red and standard elements blue. 


Reynolds-shear stresses m'm, are heavily over-predicted for the cases PhC and PhFNWM at the crest of the hill and 
in the shear layer between the recirculation region and the bulk flow. Refinement leads to an excellent match with 
reference data for the case PhF. Finally, the TKE distributions are only predicted accurately everywhere with PliF 
while PhC and PhFNWM over-predict its magnitude inside the recirculation bubble. 

The excellent results observed for Ren = 10,595 motivate an application of the wall model to a higher Reynolds 
number. For this second assessment we choose the next larger Reynolds number for which reference data is available, 
which is Re H = 19,000, and consider the same discretizations labeled PhC 19, PhF 19 and PhFNWM 19, respectively. 
The results in figure [l2| include mean velocities it\ and U 2 as well as Reynolds shear stresses «' u' 2 and are compared 
to experimental data Rapp EXP19, which do not include TKE statistics. The quality of the fine mesh including wall 
modeling is very similar to the Reynolds number Ren = 10,595 discussed above and excellent throughout and also 
matches the reattachment length perfectly. The coarser mesh PhC shows slightly worse predictions in the recircu¬ 
lation, reattachment and recovery region for the mean velocity u\ and in the recirculation bubble for the Reynolds 
stresses. Also the reattachment length is predicted significantly too short which may be due to an overly coarse mesh. 
The case without wall model PhFNWM 19 exhibits results of quality between the coarse and fine wall-modeled simu¬ 
lations. Here, another defect is highlighted: significant oscillations in the mean velocity profiles are visible especially 
at the lower wall and in the vicinity of the hill. Such oscillations are not visible for the wall-modeled computations 
and show another advantage of our wall model. 

From the investigations of flow over periodic hills the conclusion may be drawn that the present enrichment- 
based wall model exhibits favorable characteristics with respect to separated flows as well as under adverse pressure 
gradients. In this flow configuration, very coarse meshes may be used, resembling the observations made for turbulent 
channel flow. 

5.3. Backward-facing step flow 

We assess our wall modeling approach further with flow over a backward-facing step at Re /, = 5,000 with an 
expansion ratio of ER = 1.2 as studied experimentally by Jovic and Driver |62|| . DNS data of a similar configuration at 
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Figure 14: Skin-friction (left) and pressure (right) coefficients for the flow over a backward-facing step at Rei, = 5,000. 


a Reynolds number of Reh = 5,100 has been presented by Le, Moin and Kim (63j and in the context of wall modeling 
results have been presented e.g. by Chen et al. 0 mentioned earlier, who encountered difficulties predicting the 
correct skin friction and reattachment point for this flow as well. 

The computational domain behind the step is of the dimensions 30/; x 6/; x 3/; in streamwise, wall-normal and 
spanwise direction, respectively. The domain extends 30/; upstream of the step and the velocity is prescribed at 
the inflow boundary using mean DNS data of a turbulent boundary layer at a similar Reynolds number lf64l with 
an additional random perturbation of 10% of the center line velocity u c . The inflow data is only prescribed on the 
standard space u h while the enrichment u h is set to zero for simplicity. Periodic boundary conditions are applied in 
spanwise direction and slip boundary conditions at the top. The domain is meshed uniformly with four elements per 
step height /; in all space dimensions and for the case with wall modeling, three rows of elements at the lower wall 
are enriched, including the inflow region as well as the step. The resulting mesh is extremely coarse and displayed in 
figure [13] along with a contour plot of the instantaneous velocity. For the statistical results presented in the following, 
the quantities are sampled for 5,000 time steps starting after the initial transient. 

An overview over the results discussed is provided in table [3] The computation including wall modeling is labeled 
as BFS .WM 3 and compared to the same mesh where the enriched elements are replaced by standard elements, labeled 
BFS NWM. We compare our results with the experiments by Jovic and Driver ll62l labeled as BFS J&D EX/ 1 and 
the skin friction is in addition evaluated against the DNS data by Le, Moin and Kim l63ll BFS _ LMKJJNS. 

Again we start with a discussion of the distribution of the skin friction as well as the pressure coefficient along 
the lower wall. They are defined similar to d50| and (51 1 with the reference velocity u c and the reference pressure 
located at x\ = 24/;. From the results in figu re| 1 4] it can be seen that the skin friction matches reference data very well 
for the case BFS .WM3 with wall modeling. However, the computation without wall model BFS JVWM does not 
give physically reasonable results. The peak in negative skin friction is very large and shifted in streamwise direction 
by several step heights and, additionally, significant oscillations are observed. Accordingly, the reattachment length 
defined as the zero-crossing of the skin friction coefficient is predicted as X\ reatt — 6.78/; for the case with wall 
modeling which matches the references of x\j eatt = 6.0/; and x\j eatt = 6.28/; quite well. In contrast, the simulation 
without wall modeling predicts x\ jea „ = 13.49/;. An overview with respect to reattachment lengths for the simulations 
and reference data considered is also given in table [3] The prediction of the pressure coefficient shows similar quality 
as for the friction coefficient. Including wall modeling, the curve follows reference data closely, while the one without 
wall model is delayed by several step lengths. 

The mean streamwise velocity, root-mean square velocity fluctuations of the streamwise and wall-normal compo¬ 
nents and Reynolds shear stresses are displayed at six stations in figure [15] From these graphs it can be seen that both 
the velocity and fluctuations in front of the step are in good agreement with reference data implying that the simple 
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Figure 15: Mean velocity u\ in x\ direction, RMS values of the fluctuations and u' 2 in ^-direction as well as Reynolds shear stresses u'^u' 2 for 
the backward-facing step flow. 
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procedure of applying inflow data gives good results. For the case without wall modeling, the fluctuations are not 
reproduced correctly, however, which is probably due to the extreme under-resolution. 

The mean velocity behind the step matches reference data very well for the case BFS-WM 3 including wall 
modeling. It is mentioned here that the velocity is only post-processed on the element nodes which are connected 
with straight lines in the graph for simplicity. Therefore, the detailed velocity distribution at the second station inside 
the recirculation region is not shown in the graph. Without enrichment the result is barely physical and the size of the 
recirculation is significantly over-predicted. 

The RMS quantity u\ is predicted very well for the case with wall modeling. Even without wall modeling the 
match is quite good in the recovery region but lower quantities are observed inside the recirculation. The root-mean 
square of the fluctuation in ^-direction, u' v is also predicted well employing the enrichment wall modeling approach. 
In the recovery region, small discrepancies are visible, however. We assume here that this behavior is due to the 
coarseness of the mesh in the shear layer above the recirculation and is not directly related to the wall model. The 
reference without wall model yields insufficient predictions which are faulty already at the first station. The Reynolds 
shear stresses are predicted with acceptable accuracy at the first stations but an over estimation is observable around 
the fourth station. Without wall model, the Reynolds shear stresses are neither predicted accurately at the inflow nor 
behind the step. 

From the backward-facing step flow investigated with and without wall model in this section we find further 
evidence that our wall modeling approach gives excellent results in separated flow regimes. The method is robust 
with respect to jumps in bounding surfaces or ambiguous wall-normal vectors. Its strengths are accurate predictions 
for the skin-friction and pressure coefficients as well as mean velocity profiles with very coarse meshes, but even 
turbulence quantities are estimated well. 

6. Conclusion 

A new approach to wall modeling for turbulent incompressible flows at moderate and high Reynolds numbers has 
been proposed. It is suggested to enrich the function space of the computational method with a minor modification 
of Spalding’s law-of-the-wall such that the mean boundary layer gradient can be resolved with very coarse meshes. 
It is the nature of the numerical method applied, the extended finite element method, to select the most appropriate 
function among all functions available in its function space, in this case standard linear Lagrangian shape functions 
and the enrichment. Hence, Spalding’s law is not prescribed but offered as alternative solution in a consistent way. The 
general framework may be used for all kinds of enrichment functions and has been used in a variety of applications. 
As the enrichment represents near-wall turbulent fluctuations in an averaged sense and large eddies in the bulk of the 
flow are resolved, the method suggested may be interpreted as a detached-eddy simulation. In this respect, the authors 
are not seeking to point out the limitations of DES but to pioneer a new class of wall models likely to engage many 
ideas of previous wall modeling approaches in future developments. 

The method has been validated with the two most important flow regimes for wall-bounded turbulent flow, that 
is an attached boundary layer represented by turbulent channel flow and separated flow featuring a strong adverse 
pressure gradient present in flow past periodic hills and a backward-facing step. Turbulent channel flow has been 
evaluated for several Reynolds numbers giving rise to the following conclusions: The method enables the use of 
extremely coarse meshes where barely the largest eddies are resolved while high quality of the mean velocity is 
guaranteed. The results for the mean velocity profile have also shown to be essentially independent of the mesh 
employed. Root-mean-square values of the velocity fluctuations require slightly finer resolutions, however, which is 
due to the fact that the function space is ’’tuned” to represent the mean velocity and not the fluctuations. Flow past 
periodic constrictions and backward-facing-step flow exhibits the large potential of the presented method for many 
practical applications with high pressure gradient and under separated flow conditions. While standard wall models 
fail to predict the correct wall shear stresses, the present enrichment-based wall model predicts the skin-friction and 
pressure coefficients including separation and reattachment points accurately even with very coarse meshes. 

Acknowledgments 

Computational resources provided by the Leibniz Supercomputing Centre under the project pr83te are gratefully 
acknowledged. The authors thank Ursula Rasthofer for helpful and inspiring discussions on the turbulence modeling 


22 



approach and Stephan Jager for preparations towards the periodic hill benchmark. 


Appendix A. Derivatives of the enrichment in cartesian coordinates 

The first and second derivatives of the enrichment with respect to cartesian coordinates, which are required for 
evaluation of the Galerkin formulation (461, are obtained by applying the chain rule iteratively starting from equation 
(111. The equations are split into three groups: (i.) expressions in cartesian coordinates, (ii.) transformation to the 
wall coordinate y + and (iii.) derivative with respect to y + . As usual, indices define space dimensions i, j € { 1,2,3). 

i. The first derivative is 
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The ramp function is defined node-wise and interpolated with the standard FE expansion, allowing for straight¬ 
forward computation of its derivatives. 

ii. As the enrichment function is defined via the wall coordinate, its derivatives are transformed to y + yielding 

di//(x,t) d\jj dy + 
dy + dxi 

with 
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where and are obtained in a straight-forward manner via the standard FE expansion (15 i and ( jl7| . 
Applying the chain rule successively, the second derivative becomes 
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iii. The derivatives of i/j(x, t ) with respect to y + may be obtained explicitly with given <// as 
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